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Abstract 


The scattering properties of arbitrarily shaped microstrip patch 
antennas are examined. The electric field integral equation for a 
current element on a grounded dielectric slab is developed for a 
rectangular geometry based on Galerkin’s technique with subdomain 
rooftop basis functions. A shape function is introduced that allows a 
rectangular grid approximation to the arbitrarily shaped patch . The 
incident field on the patch is expressed as a function of incidence 
angle 0 l i <j> 1 . The resulting system of equations is then solved for 
the unknown current modes on the patch , and the electromagnetic 
scattering is calculated for a given angle. Comparisons are made 
with other calculated results as well as with measurements. 

Introduction 

In the early 1980’s, the moment method technique was developed for analyzing microstrip 
patch antennas by using the spectral domain Green’s function. This technique accurately 
accounts for dielectric thickness, dielectric losses, and surface wave losses and can be extended 
to include the effects of a cover layer with a different dielectric constant on top of the antenna. 
Because of its spectral nature, the technique can be easily extended to model an infinite array 
of patches by the examination of only a single unit cell. 

Early papers on this subject include those by Bailey and Deshpande (refs. 1 to 3) and 
Pozar (ref. 4). In those papers the authors used both subdomain and entire-domain basis 
functions to model the current on the patch. Results such as bandwidth, input impedance, 
and resonant frequency were presented for rectangular patches. The use of subdomain basis 
functions yields more flexibility in the modeling of the patch current, whereas the use of entire- 
domain basis functions yields a smaller number of unknowns in the solution. For this reason, 
many subsequent analyses involve entire-domain basis functions that are limited to canonical 
shapes such as rectangles, circles, and ellipses. Using entire-domain basis functions, Pozar and 
Schaubert (ref. 5) extended the method from a single patch to an infinite array of patches. 
Aberle and Pozar (refs. 6 and 7) analyzed circular microstrip patch antennas, in both single 
and array geometries, using entire-domain basis functions. Bailey and Deshpande (ref. 8) also 
used entire-domain basis functions in their study of elliptical and circular patches. Aberle and 
Pozar (refs. 9 and 10) have improved on the original idealized probe feed model by including 
attachment modes of current that accurately model the current singularity at the probe feed 
point. 

Recently, much work has been published regarding the scattering properties of microstrip 
antennas on various types of substrate geometries. Virtually all this work has been done with 
entire-domain basis functions for the current on the patch. Newman and Forrai (ref. 11) have 
analyzed the electromagnetic scattering from a rectangular patch. Jackson (ref. 12) extended 
this research to include a superstrate covering the antenna. Pozar (ref. 13) examined a patch 
on a uniaxial substrate. Aberle, Pozar, and Birtcher (ref. 14) included the improved probe feed 
model for analyzing the scattering from circular patches. 

Some work has been published concerning the use of subdomain basis functions for modeling 
the current on the patch antenna. Most of this work was done in the spatial domain and 
cannot be extended to infinite arrays as in the spectral domain approach. Hall and Mosig 
(refs. 15 and 16) have analyzed rectangular microstrip antennas using a mixed-potential integral 
equation with subdomain rooftop basis functions. Mosig (ref. 17) has also used this approach 



to analyze patch antennas of arbitrary shape. Michalski and Zheng (ref. 18) have used a similar 
formulation with triangular-surface patch basis functions, which are similar to finite element 
techniques, to model patches of arbitrary shape. Martinson and Kuester (ref. 19) have used a 
network approach along the edge of the patch to examine different shapes. As with the spatial 
domain method, this approach is not easily extended to array geometries. Hansen and Janhsen 
(ref. 20) have outlined a spectral domain approach that uses subdomain basis functions for 
modeling rectangular patches with a microstrip line-feed network. 

This paper describes spectral domain analyses of arbitrarily shaped microstrip patch antennas 
in which subdomain basis functions are used to model the patch current. To simplify the analyses, 
the antenna feed will not be considered. The antenna is considered to be open circuited from the 
feed network (i.e., the feed impedance is infinite). Results are presented in the form of scattering 
as a function of frequency for a few representative shapes. Comparisons are made with measured 
data and with results from other analysis techniques. 
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dyadic Green’s function 

component of the spatial domain Green’s function 
component of the spectral domain Green’s function 
amplitude of mode mn 

surface current on the microstrip patch antenna 

= V—i 

variables of integration in cylindrical coordinates 

propagation constant for free space, 27r/A 0 

spectral domain transformation variable for the x-direction 

spectral domain transformation variable for the ^-direction 

propagation constant for the dielectric slab in the z-direction 

propagation constant for free space in the z-direction 

dimension of the patch in the x-direction 

dimension of the patch in the ^-direction 

number of subdivisions in the x-direction 

number of subdivisions in the y-direction 
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Abbreviation: 

dBsm 


characteristic equation for transverse electric modes 

characteristic equation for transverse magnetic modes 

component of excitation voltage vector 

coordinates of the field point 

coordinates of the current mode rrm 

coordinates of the source point 

unit vector in the x-direction 

unit vector in the y-direction 

impedance matrix to be solved 

impedance of free space, 377 ohms 

unit vector in the ^-direction 

= v/K| + 

cell size in the x-direction 

cell size in the y-direction 

relative permittivity of the dielectric slab 

incident angle of electromagnetic wave 

piecewise linear function for the current on the patch 

permeability of free space 

pulse function for the current on the patch 

0 polarized backscatter from 9 polarized incident field 

frequency of the electromagnetic field, rad 

a unit denoting decibels referenced to square meters 


Theory 


The geometry of a rectangular microstrip patch antenna is shown in figure 1. The patch is 
on a grounded dielectric slab of infinite extent. The dielectric slab has a relative permittivity e r 
and thickness d. Assuming that the patch is perfectly conducting, the boundary condition on 
the patch is given by 

E inc tp scat /i ^ 

tan — “Man 

The incident field is the field at the patch location attributable either to an incident plane wave 
or a probe or stripline feed. The scattered field is found from the currents excited on the patch 


as 


gscat 


(x,y,z) 



, y, z\x 0 , y 0 , z 0 ) ■ J {x 0 , y Q , z 0 ) dx 0 dy 0 dz 0 


( 2 ) 


<-> 

where G is the dyadic Green’s function for a current element on a grounded dielectric slab and 
J is the electric current density for the unknown vector on the patch. The dyadic Green’s 
function can be written as 


G - xG XI x + xGxy y + xG xz z + yGy X x + yG yy y + yG yz z + z G zx Si + zG zy y + zG zz z (3) 
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where 


I roc roc _ 

G “ 6 = 4 tt 2 J j G(lb ^ Kx ' Ky ' z ^ z °) exp ^ Kx ^ x ~ x ° ^ exp ^ Ky ( y ~ Vo ^ dKxdK y ( 4 ) 


and a and b can be x, y , or 2 . Note that the components for the Green’s function G a b are known 
in the spectral domain and must be transformed back to the x, y domain; hence, the two infinite 
integrals in equation (4) are required. 

The components of the Green’s function are given by 
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where 


T m = e r K 2 COs(Kid) + jKi sin (K\d) 


(ii) 

T t = K\ cos (Kid) + j K‘i sin ( K j d) 


(12) 

K\ = yJe r Kl - 0* 

Im (K\) < 0 

(13) 

k-2 = \JkI- [K 

Im (K 2 ) < 0 
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* 

+ 
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The remaining terms of the Green’s function are not needed in this analysis. Details of the 
derivation of the Green’s function can be found in reference 3. Additional forms of the Green’s 
function that include the effects of a dielectric cover layer above the antenna are available in 
reference 10. 
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The current density J is modeled as a summation of piecewise linear subdomain basis 
functions known as rooftop basis functions. This approach contrasts with use of the entire- 
domain basis functions that span the entire patch. Entire-domain basis functions, such as sines 
and cosines, are useful for analyzing rectangular or circular patches, but become cumbersome for 
other shapes. Mathematically, the subdomain basis functions for the components of the current 
are described as 

M N -\- 1 

j,(*, i)=ee ir A™ <x) n„ (») (is) 

m = 1 n— 1 


Af+1 N 

J y (x, y) = Y, E r r A n (v) n m (x) 
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where the functions A and II are “triangle” and “pulse” and are expressed as 
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1 
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Otherwise J 


(17) 


(18) 


(19) 


where Ax = 2L X /{M + 1) and Ay = 2 L y /(N + 1). 

When equations (2) and (4) are combined, the order of integration may be changed and the 
basis functions that represent the patch current density may be taken into the transform domain. 
These current density functions for the spectral domain are given by 

M N + 1 
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where 


F™ n (K X1 K y ) = Ax Ay 


sin (KyAy/2) 
KyAy/2 
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K x Ax/2 
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F™ (K x , K y ) = Ax Ay 


sin ( KyAy/2 ) 
KyAy/2 


sin (/S'* Ax/2) 
K/Ax/2 


exp [ j K x x m jK y y n + jK x (Ax/2)] 


(23) 


Galerkin’s method can be applied to the resulting equations to test them with the same set 
of basis functions. The solution yields a set of simultaneous equations that can be solved with 
standard techniques. Symbolically, this approach is represented as 


JJ J pq ■ E l t "n dxdy = -JJ/ Pq - E tan dxdy 


(24) 
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Note that the integration is over x and y instead of x Q and y Q . On the right side of equation (24) 
the x, y integration may be performed. The resulting Fourier transforms are similar to those 
described in equations (22) and (23). These equations can be shown in matrix notation as 
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where the impedance matrix terms are given by 

j roc roo ~ 
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The integrations in equations (26) (29) must be done numerically but can be simplified with 
the following change of variables: 


K x = K cos a K y = K sin a (30) 

With this change of variables, the integrals are changed to the form 



] dK x dK y 



] Kd Kda 


(31) 


The integration from 0 to 27T may be further reduced to an integration from 0 to 7r/2 based on 
the even and odd properties of the integrand. Each of the four submatrices in the impedance 
matrix is of Toeplitz form, so only the first row of each submatrix needs to be calculated by 
numerical integration. The remaining terms can be filled in with these terms. Furthermore, 
because the impedance matrix terms Z§y mTl = . even more computer time is saved. 

To examine the scattering from a microstrip patch antenna, the left side of equation (25) 
must be evaluated. Each member of the excitation vector can be written as 

V Pq = Jis JW ' EmC dxdy (32) 

which is the incident field reacted with each pq current mode on the patch. After we use 
reciprocity, equation (32) can be rewritten as 


vPq _ ■ E ° 
jupo 


( 33 ) 
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In equation (33), E 0 is the vector amplitude of the incident plane wave, W q is the far- field 
radiation from vector current mode pq on the patch, and /jup 0 is the required strength 
of an infinitesimal dipole source to produce a unit amplitude plane wave. The incident plane 
wave is from the direction 0 ? , <f) 1 in spherical coordinates with components Eq and E ( p. Typical 
scattering results are of the form 


a ee = 4nr 2 |£| cat | 2 ( 34 ) 

which is the 0 polarized backscatter from a unit amplitude 0 polarized incident field. 

The fields radiated by a current mode on the patch can be found from the Green’s function 
(see eq. (5)). The field at the point x , y, z from an x directed source located at the point x 0y y 0 , d 
is given by 


Eb(x,y,z) 



exp \jK x (x 


x 0 )] exp \jK y (y — y 0 )\ exp [- jK 2 (z - d)] dK x dK y 


(35) 


where b can be either x, y, or z. Likewise, the values x,y,z from a y directed source at point 
x 0 ,y 0 ,d are given by 


1 yoo r x __ 

E b {x,y,z )= — -2 / / Ghy exp[jK I (x - Xo)] cxp[jKy(y -y a )] exp[-jK 2 (z - d)]dK x dK y (36) 

J — 00 j -DC 

and again b can be either x, y, or z. These equations can be evaluated by the method of 
stationary phase and integrated over the extent of each basis function to give the fields radiated 
by that basis function in the presence of the grounded dielectric slab (ref. 13). After this 
evaluation has been done and the resulting equations are converted to spherical coordinates, the 
far-field components due to a single x directed current mode are 


E T M,0 = 


Zo 

271 - 


exp (~jK 0 r) 


( . js j\ q E\ E 0 cos (j) sin (Aid) ^77^ ^ ^ \ 

exp (jK 2 d) cos 0 — F x (A*, K y ) 
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Er^9:4>)=^ 


2i r 


exp (—jK 0 r) 


, n n —Kn sin d> sin (Aid) ,,, _ _ x 

exp (jK 2 d) cos 9 — F x ( K x , K y ) 


(38) 


where K x and K y are evaluated at the stationary phase points 

K x = — K 0 sin 9 cos <p 
Ky = — K 0 sin 6 sin 0 
Similarly, the fields radiated by a single y directed current mode are given by 
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(39) 
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exp (jK 2 d) cos 6 


cos <f> sin (K\d) 
T e 


F? n {Kx ,Ky ) 
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where K x and K y are the same as in equation (39). By using equations (37) to (41) in 
equation (33), we can determine the left side of equation (25). 

After the impedance matrix and the excitation vector have been calculated, the simultaneous 
equations can be solved for the unknown current coefficients. Then, the scattered fields can be 
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calculated by a summation of the radiated fields from each mode on the patch. If the patch is 
rectangular, this process is straightforward. However, if the patch is some other shape, additional 
steps are needed to model it properly. 

Consider the irregular patch shown in figure 1. To predict the scattering from this patch, first 
enclose it within a rectangle. The impedance matrix and excitation vector can be calculated for 
this rectangular patch. A shape function is introduced that is equal to 1 for each mode that 
has its center point inside the irregularly shaped patch and equal to 0 if the center point of 
the mode is outside the irregularly shaped patch. The set of simultaneous equations can also 
be modified to consider only the modes for which the shape function is equal to 1. Thus, the 
boundary of the irregularly shaped patch is approximated by a rectangular grid or “stair step.” 
As the number of subdivisions increases, the approximation to the true boundary of the patch 
improves. However, as the number of subdivisions increases, the time required to compute and 
solve the impedance matrix increases as well. 

Results 

Computer programs have been written to solve the matrix equation (25). These programs 
are listed in the appendix. The first program computes the elements of the impedance matrix 
by numerical integration. As mentioned previously, only the first row of each submatrix is 
calculated. The rest can be filled in by rearranging the first row. Also, because the Z xy and 
Zy X submatrices are related, only the Z X y submatrix is calculated. The impedance matrix is 
then stored in a data file. The second program reads in the impedance matrix from the file and 
calculates the excitation vector for the given angle of incidence. The system of equations is solved 
and the electromagnetic scattering is calculated at the same angle. If scattering information is 
required over a band of frequencies, the third computer program reads in and arranges impedance 
matrices for several frequencies and computes the scattering at a given incidence angle as a 
function of frequency. It is necessary to compute only the impedance matrix at a few widely 
spaced frequencies because the impedance matrix terms are slow to change as a function of 
frequency. The impedance matrix for other frequencies can be found by an interpolation of each 
impedance matrix element. Newman and Forrai (ref. 11) have used this approach successfully 
with entire-domain basis functions. 

To ensure that the computer programs are correct, comparisons are shown in figure 2 for 
the calculated and measured data presented by Newman and Forrai (ref. 11) and the calculated 
results from the subdomain technique. Because both the entire-domain basis functions used in 
reference 1 1 and the subdomain basis functions described here model the patch shape accurately, 
the calculated results from each technique should agree. The number of subdivisions here was 
chosen to be M = N — 12, which is adequate for modeling the patch across the frequency 
band. If frequency were further increased, more patch subdivisions would be necessary. Figure 2 
shows the impedance matrix that was calculated through numerical integration at frequency 
steps of 400 MHz and that was interpolated at frequencies between these steps. A quadratic 
interpolation technique was used here as well as for the following results. Close agreement 
between the calculated results can be seen in figure 2; the only slight differences are at the peaks 
of the curves. Also, a slight offset in frequency is noted between the two calculations. However, 
this offset is not uncommon when two techniques are compared and is likely attributable to 
minor differences in the computer codes. The measured results shown in figure 2 agree well 
with the calculated results in some areas and disagree in others. Again, a slight frequency shift 
is noted when compared with the calculated results; this shift may indicate physical tolerances 
of the patch size, substrate thickness, or substrate dielectric parameters. At some points, the 
measured data agree better with the data from the subdomain technique; at other points, the 
measured data agree better with those from the entire-domain technique. In some areas, the 
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measured data do not agree with any calculated results; rapid fluctuations of the measured data 
in these instances suggest a problem with the measurement or calibration of the radar data. 

Next, the shape function was included in the computer programs and two different circular 
microstrip patch antennas were modeled. In figure 3 calculated results for a circular patch 
modeled with subdomain basis functions are compared with results calculated by Aberle and 
Pozar (ref. 6) with entire-domain basis functions. The impedance matrix case for the subdomain 
basis function was calculated at 500-MHz steps and was interpolated elsewhere. For entire- 
domain basis functions, the patch boundary is a perfect circle; for the subdomain basis functions, 
the patch boundary only approximates a circle. As before, a slight frequency shift is noted 
between the two approaches. This shift can be attributed to slight differences between the two 
computer codes and the subdomain approximation of a circular boundary. The frequency shift 
increases as the frequency increases and is expected, as the difference between the true and 
approximated boundary is larger electrically as frequency increases. If the frequency shift is 
ignored, the two results are similar across the frequency band, especially from 2 to 4.5 GHz. 
Approximation of a circular boundary in the subdomain calculation may explain the slight 
differences in the relative power levels. 

A second circular patch was modeled and compared with results from entire-domain calcu- 
lations and with measured data. The entire-domain calculations and the measured data were 
supplied by Aberle, Pozar, and Birtcher (ref. 14) and are shown in figure 4. The agreement 
is fairly good between the subdomain, entire-domain, and measured data. Again, a slight fre- 
quency shift in the data is noted, as are slight differences in the values at the resonant peaks. 
For both patches, the number of subdivisions was chosen to be M = N = 12. More subdivisions 
were used for the latter patch to improve the agreement, but little improvement was noted. 

Two other shapes, an equilateral triangle and a trapezoid, have been modeled with subdomain 
basis functions. As before, the impedance matrices were calculated at 500-MHz steps and were 
interpolated at other frequencies. For these two shapes no entire-domain basis function results 
have been reported. The measured data were collected in the Experimental Test Range at the 
Langley Research Center. 

The calculated and measured results for the triangular microstrip patch antenna are shown in 
figure 5. As before, the number of subdivisions was chosen to be M — N = 12. As in all previous 
cases, a slight frequency shift is evident in the data and the peak values of the data are slightly 
different. The same rapid fluctuations in the measured data as in figure 2 are observed from 
9 to 13 GHz; these fluctuations are caused by imperfections in the measurement and calibration 
process. The measured versus calculated results for the trapezoidal microstrip patch are shown 
in figure 6. Because the trapezoid is larger than the earlier examples, the number of subdivisions 
of the patch was changed to M = 10 and N = 20. Other combinations of M and N were tried 
with little change in the results. 

The results in this case are not as good as in the previous cases. The shift in frequency 
between the calculated and measured data shown in figure 6 is larger than those in the previous 
cases. The relative power levels, however, are quite close for most peaks, although a large 
difference is evident at approximately 7.3 GHz. Rapid fluctuations in the measured data at that 
point suggest a measurement problem as the cause of the discrepancy. 

Conclusions 

A subdomain moment method technique has been developed to examine the scattering prop- 
erties of microstrip patch antennas. Antennas of virtually any shape may now be analyzed with 
this technique merely by changing the shape function that defines the outer antenna boundary. 
Results were presented for rectangular, circular, triangular, and trapezoidal microstrip patch an- 
tennas. In all cases, slight shifts in frequency between subdomain, entire-domain, and measured 
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data were noted, as were slight differences in power level at the response peaks. The sub- 
domain calculations have, however, predicted the correct general scattering from all the patches 
examined. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
September 30, 1992 
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Appendix 

Computer Programs That Solve Matrix Equation (25) 

Evaluation of Z xx , Z yy , and Z xy Submatrices 

The following program computes the first row of the Z xx , Z yy , and Z xy submatrices by 
numerical integration. The remaining terms of the impedance matrix are then filled and the 
matrix is stored in a data file. 

ccccc Patch Impedance Matrix 
ccccc Subdomain Basis Functions 
ccccc 

program IMPEDANCE 

parameter (mm=12 ,nn=12,mmpl = (mm+l) ,nnpl=(nn+l) ) 

parameter (nmaxl= (mm*nnpl) ,nmax2=(mmpl*nn) ,nmax=(nmaxl+nmax2) ) 

real x(mmpl) ,y(nnpl) , Lx,Ly ,dx ,dy ,dx2 ,dy2,Dxy ,d 

real Ko,fr,pi,pi2,sl ,s2,epo,muo 

real B2 , A2 , DELTY2 , FACTY2 , FACTY1 , FY2 

real alpha, ca, sa, deltyr, deltyi , YI1 , YI2 , ANGC60) 

integer NQ2 ,NUM,NS (60) ,Con(3) 

complex cj , er ,Ke , tempx , tempy ,P(60) ,FY1 

complex Gxx , Gxy , Gyy , sxx , sxy , syy , Txx , Tyy , Txy 

complex K,K1 ,K2 ,Tm Je ,N1 ,N2 ,D1 ,D2 ,C1 , sd 

complex ZZ(nmax ,nm ax) , DELTYI 

complex zxx(mmpl ,nnpl ,mmpl ,nnpl) , zxy (mmpl ,nnpl ,mmpl ,nnpl) 
complex zyx(mmpl ,nnpl ,mmpl ,nnpl) ,zyy(mmpl ,nnpl ,mmpl ,nnpl) 
complex axx(mmpl ,nnpl) ,axy (mmpl ,nnpl) ,ayy(mmpl ,nnpl) 
complex bxx(mmpl ,nnpl) ,bxy (mmpl ,nnpl) ,byy (mmpl ,nnpl) 

DIMENSION U1 (3) ,U2(10) ,R1(3) ,R2(10) ,U(13) ,R(13) 

EQUIVALENCE (Ul(l) ,U(1)) , (U2(l) ,U(4)) , (Rl(l) ,R(1)) , (R2(l) ,R(4)) 

DATA Ul/. 11270166537925, .5, . 88729833462074/ ,U2/ . 01304673574141 , .06 

1746831665550. . 16029521585048. . 28330230293537 42556283050918 . .5744 

23716949081. . 71669769706462. . 83970478414951 . . 93253168334449 . .986953 
326425858/ ,R1/ . 27777777777777 , . 44444444444444 , . 27777777777777/ ,R2/ . 

403333567215434. . 07472567457529. . 10954318125799 .. 13463335965499 . .14 

5776211235737. . 14776211235737. . 13463335965499 10954318125799 . .0747 

62567457529. . 03333567215434/ 

CCCCC 

ccccc ALL DIMENSIONS IN METERS 
ccccc 

cj=(0. 0,1.0) 

pi=2 . 0*asin(l . 0) 

pi2= pi/2.0 

epo= 8 . 854E-12 

muo= pi* (4 . 0E-07 ) 

er= cmplx(2.2,-0.0022) 

fr= 6.0E+09 

Ko= 2. 0*pi*fr*sqrt (muo*epo) 

Ke= 2.0*pi*fr*csqrt (muo*epo*er) 
d= 0.00159 
Lx= 0.023 
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Ly= 0.023 

dx= (2 . 0*Lx) / (float (mmpl) ) 

dy=(2.0*Ly)/(f loat (nnpl)) 

dx2=dx/2 . 0 

dy2=dy/2.0 

Dxy= dx*dx*dy*dy 

do 3 i=l,mmpl 

x(i) = -Lx + (f loat (i) *dx) 

3 continue 
do 5 i=l,nnpl 
y(i) = -Ly + (f loat (i) *dy) 

5 continue 

write (6 , 10) Ko ,Ke ,Lx ,Ly 

10format(’Ko=’ ,f7.3, ’ Ke= * ,f7 . 3 ,f 7 . 3 , ’ Lx=’,f9.6,’ Ly=’,f9.6) 
zxx=(0 .0,0.0) 
zxy=(0 .0,0.0) 
zyx=(0 .0,0.0) 
zyy= (0 . 0 , 0 . 0) 
axx= (0 . 0 , 0 . 0) 
axy=(0 ,0,0.0) 
ayy=(0 .0,0.0) 
ccccc 

ccccc Calculate Integrals for the Z matrix 
ccccc 

ccccc Define Limits for K integration 
ccccc 

P(l)= (0.0,0. 0) 

P(2)= cmplx( (0 . l*Ko) , (0 . l*Ko) ) 

P(3)= cmplx((1.0*Ko),(0.1*Ko)) 

P(4)= cmplx( (1 . l*real (Ke) ) , (0 . l*Ko)) 

P(5)= cmplx((l . l*real (Ke) ) , (0 . 0) ) 

P (6) = cmplx ( ( 1 . 5*real (Ke) ) ,0.0) 

P(7)= cmplx( (4. 0*real (Ke) ) ,0.0) 

P(8)= cmplx((5,0*real(Ke)) ,0.0) 
do 11 n=9,60 
P(n)= f loat (n-8) *6000.0 
11 continue 
NS=50 
NS(1)=10 
NS (2)=20 
NS (3) =20 
NS(4)=10 
NS (5) =10 
NS (6) =40 
NS(7)=20 
NS (8)=100 
ccccc 

ccccc THETA Integration first 
ccccc 

B2=pi2 
A2=0 . 0 
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NQ2=6 

DELTY2= (B2-A2) /FLOAT (NQ2) 

C+FIND INITIAL VALUE FDR THETA INTERVAL. 

DO 200 11=1, NQ2 
YI2 = float (II - 1) 

FY2 = A2 + (YI2*DELTY2) 

CEVALUATE FUNCTIONS AT 10 POINTS PER THETA INTERVAL. 

DO 180 L2=l , 10 
FACTY2=R (3+L2) 

alpha=FY2 + (U (3+L2) *DELTY2) 
ca=cos (alpha) 
sa=sin (alpha) 

C 

C NOW START K INTEGRATION 
C 

Con=0 

Txx=(0. 0,0.0) 

Txy=(0. 0,0.0) 

Tyy=(0. 0,0.0) 

DO 150 IS=1 , 59 
bxx=(0 ,0,0.0) 
bxy=(0 ,0,0.0) 
byy=(0.0,0.0) 

DELTY1= (P(IS+1) -P ( IS) ) /FLOAT (NS (IS) ) 

C*FIND INITIAL VALUE FOR K INTERVAL. 

DO 100 JJ=1 ,NS(IS) 

YI1 = f loat ( JJ - 1) 

FY1 = PUS) + (YI1+DELTY1) 

CEVALUATE FUNCTIONS AT 10 POINTS PER K INTERVAL. 

DO 80 Ll=l , 10 
FACTY1=R(3+L1) 

C 

C Find K and evaluate Green's functions 
C 

K=FY1 + (U (3+L1) *DELTY1) 

Kl=csqrt((Ke**2)-(K**2)) 
if (aimag(Kl) .gt.O.O) K1 = conjg(Kl) 

K2=csqrt ( (Ko**2) - (K**2) ) 

if (aimag(K2) .gt.O.O) K2 = conjg(K2) 

sd=csin(Kl*d) 

Te= (Kl*ccos(Kl*d) ) + (cj *K2*csin(Kl*d) ) 

Tm= (er*K2*ccos(Kl*d)) + (cj *Kl*csin(Kl*d) ) 

Cl= ( (cj *377 . 0*sd)/ (Ko*Te*Tm) ) 

Gxx= (-C1) * ( (Kl*K2*ca*ca*Te) + (Ko*Ko*sa*sa*Tm) ) 
Gyy= (-C1) * ( (Kl*K2*sa*sa*Te) + (Ko*Ko*ca*ca*Tm) ) 
Gxy= Cl * ( (Ko*Ko*ca*sa*Tm) - (Kl*K2*ca*sa*Te) ) 

C 

C Evaluate Basis functions for each mode 
C 

Dl= dx2*K*ca 
Nl= csin(Dl)/Dl 
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D2= dy2*K*sa 
N2= csin(D2)/D2 
C 

C multiply by K for polar integration 
C 

sxx= Gxx*K*FACTY 1 *Dxy* (Nl**4) * (N2**2) 
sxy= Gxy*K*FACTYl*Dxy* (Nl**3) * (N2**3) 
syy= Gyy*K*FACTY 1 *Dxy* (Nl**2) * (N2++4) 

C 

C SET p and q = 1 and then vary m and n 
C 

do 15 m=l ,mmpl 

tempx= ccos(K*ca*(x(m)-x(l))) 
do 12 n=l,nnpl 

tempy= ccos(K*sa* (y (n)-y (1) ) ) 
bxx(m,n)=(sxx*tempx*tempy) + bxx(m,n) 
byy (m,n)=(syy*tempx*tempy) + byy(m,n) 

12 continue 
15 continue 
C 

C Take abs() of dx and dy terms and do the sign later 
C 

do 25 m=l > mmpl 

tempx= csin(K*ca*abs(x(l)-(x(m)-dx2) )) 
do 20 n=l,nnpl 

tempy= csin(K*sa*abs (y (1) -dy2-y (n) ) ) 
bxy(m,n)= (-1 . 0*sxy*tempx*tempy) + bxy(m,n) 

20 continue 
25 continue 

ccccc write (6,75)K, (sxx/FACTYl) 

75 format (f 12.5 ,f 12.5 ,f 15. 12,f 15. 12) 

80 continue 
100 continue 

Txx=(bxx(l , 1) *DELTY1) + Txx 
Txy= (bxy (1,1) *DELTY 1 ) + Txy 
Tyy= (byy (1,1) *DELTY 1 ) + Tyy 
axx= (bxx*DELTY 1 *FACTY2) + axx 
axy= (bxy +DELTY 1 *FACTY2) + axy 
ayy= (byy*DELTYl*FACTY2) + ayy 
ccccc write (6 , 120)P(IS+1) , (bxx(l , 1) *DELTY1) ,Txx 

120 f ormat (f 10 . 2, f 10 . 2,f 13 . 10, f 13. 10, f 13. 10, f 13. 10) 

if (cabs(bxx(l , 1)*DELTY1) .It. (0 . 02*cabs (Txx) ) ) Con(l)=l 
if (cabs (bxy (1 , 1)*DELTY1) . It. (0.02* cabs (Txy))) Con(2) = l 
if (cabs (byy (1 , 1) +DELTY1) . It . (0 . 02*cabs(Tyy) ) ) Con(3)=l 
if ( (Con(l) . eq. 1) . and. (Con (2) . eq. 1) . and. (Con (3) . eq. 1) 
&,and. (IS.gt .16)) goto 170 
150 continue 
170 continue 

ccccc write (6 , 175) (alpha* 180 .0/pi) ,P(IS+1) 

175 format (f 12 . 2, f 12. 2, f 12. 2) 

180 continue 
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200 continue 

202 format (f 8.3, f 8.3, f 8.3, f 8.3, f 8.3, f 8. 3) 


C 

C Multiply by 4 for 1 quadrant integration 
C Divide by 2pi*2pi for inverse Fourier transform 
C Multiply by -1 for E in = - E scat 
C 

axx= ( (-1 . 0/(pi*pi) ) *axx*DELTY2) 
axy= ( (-1 . 0/ (pi*pi) ) *axy*DELTY2) 
ayy= ((-1.0/ (pi*pi) ) *ayy*DELTY2) 

C 

C Now arrange the other matrices 
C 

do 310 ip=l,mm 
do 310 iq=l,nnpl 
irow= ( (ip-1) *nnpl) + iq 
do 300 m=l,mm 
do 300 n=l,nnpl 
icol= ((m-l)*nnpl) + n 
if ( (ip . eq. 1) . and. (iq. eq. 1) ) then 
zxx(ip, iq,m,n)= axx(m,n) 
else 

mp=iabs (ip-m)+l 
nq=iabs (iq-n) +1 
zxx(ip , iq,m,n)= axx(mp,nq) 
endif 

ZZ(irow, icol)= zxx(ip , iq,m,n) 

300 continue 
310 continue 
C 

do 360 ip=l,mm 
do 360 iq=l,nnpl 
irow= ( (ip-1) *nnpl) + iq 
do 350 m=l,mmpl 
do 350 n=l,nn 

icol= ((m-l)*nn) + n + nmaxl 
sl=sign(l .0, (x (ip) - (x (m) -dx2) ) ) 
s2=sign(l .0, (y (iq) -dy2-y (n) ) ) 
if ( (ip. eq. 1) . and. (iq. eq. 1) ) then 
zxy(ip , iq,m,n)= axy (m,n) *sl*s2 
else 

mp=iabs (ip-m) +1 
if (m.lt.ip) mp=mp+l 
nq=iabs (iq-n)+l 
if (n . It . iq) nq=nq-l 
zxy(ip, iq,m,n)= axy (mp,nq) *sl*s2 
endif 

ZZ(irow, icol)= zxy (ip , iq,m,n) 
zyx(m,n, ip, iq) =zxy (ip , iq,m ,n) 

350 continue 
360 continue 
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c 


do 380 ip=l,mmpl 
do 380 iq=l,nn 

irow= ((ip-l)*nn) + iq + nmaxl 
do 375 m=l ,mm 
do 375 n=l,nnpl 
icol= ( (m-1) *nnpl) + n 
ZZ(irow, icol)= zyx(ip, iq,m,n) 

375 continue 
380 continue 
C 

do 400 ip=l,mmpl 
do 400 iq=l,nn 

irow= ((ip-l)*nn) + iq + nmaxl 
do 400 m= 1 , mmp 1 
do 400 n=l,nn 

icol= ((m-l)*nn) + n + nmaxl 
if ( (ip.eq. 1) . and. (iq. eq. 1) ) then 
zyy(ip. iq,m,n)= ayy(m,n) 
else 

mp=iabs (ip-m)+l 
nq=iabs (iq-n)+l 
zyy (ip, iq,m,n)= ayy(mp,nq) 
endif 

ZZ(irow, icol)= zyy (ip , iq,m,n) 

400 continue 
C 

ccccc do 1200 m= 1 , nmax 

ccccc do 1200 n=l,nmax 

ccccc write (6 , 1100)m,n, ZZ(m,n) 

ccccc 1100 f ormat (i4 , i4 , f 12 . 8 , f 12 . 8) 
ccccc 1200 continue 

open ( 12 , FILE= * IMP 9 , STATUS= J NEW > , F0RM= J UNFORMATTED ' ) 
write ( 12) fr ,Lx,dx,Ly,dy ,d,Ko,Ke ,er ,x,y 
write(12)ZZ 
close(12) 

2000 stop 
end 
ccccc 

Solution for a Given Incidence Angle 

The following program reads an impedance matrix data file, computes the excitation vector, 
and computes backscatter for a given incidence angle. The shape function is introduced and 
only those modes with a shape function equal to 1 are retained in the solution. 

ccccc Microstrip Patch Matrix Inversion Code 

parameter (mm-12,nn=12,mmpl=(mm+l) ,nnpl= (nn+1) ) 

parameter (nmaxl=(mm*nnpl) ,nmax2=(mmpl*nn) ,nmax= (nmax 1+ nmax 2) ) 

real Lx,dx ,Ly ,dy ,Ko ,fr ,x(mmpl) ,y(nnpl) ,px,py 

real pi , dx2 , dy2 , Theta , Phi , Kx , Ky , K2 , FI , F2 

integer ipvt (nmax) ,Sx(mm,nnpl) ,Sy(mmpl ,nn) 
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complex er , Ke , axx , bxx , r es , Rs , Di , Et , Ep , Ph 

complex c j , K1 , Tm , Te , Ptx , Pty , Ppx , Ppy , Etx , Ety , Epx , Epy 

complex fac(nmax,nmax) ,wk(nmax) 

complex ZZ(nmax ,nmax) ,V(nmax) ,Vx(nmaxl) ,Vy(nmax2) 
complex ZZN (nmax , nmax) , VN (nmax) , Vmid , Jmid 
complex J(nmax) 
character*8 fn 

equivalence (Vx(l) ,V(1)) , (Vy(l) ,V(nmaxl+l)) 

DIMENSION U1 (3) ,U2(10) ,R1(3) ,R2(10) ,U(13) ,R(13) 

EQUIVALENCE (U1 (1) ,U(1) ) , (U2(l) ,U(4) ) , (R1 Cl) ,R(1) ) > 0*2(1) ,R(4)) 
DATA Ul/. 11270166537925, .5, . 88729833462074/ ,U2/ .01304673574141 , .06 

1746831665550. . 16029521585048. .28330230293537. .42556283050918. .5744 

23716949081. . 71669769706462. .83970478414951. .93253168334449. .986953 
326425858/ , R1 / . 27777777777777 , . 44444444444444 , . 27777777777777 / , R2/ . 

403333567215434 . . 07472567457529 . . 10954318125799. . 13463335965499 . . 14 

5776211235737. . 14776211235737. . 13463335965499. . 10954318125799. .0747 

62567457529. . 03333567215434/ 

C 

read(5, l)fn 

1 format (A) 

open( 12 , f ile=f n , status= ’ OLD ’ , f orm= ’UNFORMATTED ’ ) 

read ( 1 2 ) f r , Lx , dx , Ly , dy , d , Ko , Ke , er , x , y 

read(12)ZZ 

close(12) 

wr ite (6 , 2) f r , Lx , Ly , d , Ko , Ke , er 

2 format (elO . 3,f 10 . 6,f 10 .6,f 10. 6,5f 10.4) 
dx2=dx/2 . 0 

dy2=dy/2 .0 

pi= 2.0*asin(l .0) 

cj= (0.0, 1.0) 

CCCCC 

CCCCC Find Vx and Vy excitation for angle Theta and Phi 
CCCCC 

Theta= 60. 1* (pi/180.0) 

Phi= 0.1*(pi/180.0) 

Et= 1.0 
Ep= 0.0 

Kx= -Ko*sin(Theta)*cos(Phi) 

Ky= -Ko*sin(Theta) *sin(Phi) 

Kl= csqrt ( (Ke**2)-(Kx**2)-(Ky**2) ) 
if (aimag(Kl) .gt.0.0) Kl= conjg(Kl) 

K2= Ko*cos (Theta) 

Tm= (er*K2*ccos (Kl*d) ) + (cj *Kl*csin(Kl*d) ) 

Te= (Kl*ccos (Kl*d) ) + (cj*K2*csin(Kl*d) ) 

Ptx= (Kl*Ko*cos (Phi) *csin(Kl*d) ) /Tm 
Pty= (Kl*Ko*sin (Phi) *csin(Kl*d) ) /Tm 
Ppx= - (Ko*Ko*sin(Phi) *csin(Kl*d) ) /Te 
Ppy= (Ko*Ko*cos(Phi) *csin(Kl*d) )/Te 
ccccc 

Di= - (4 . 0*pi) / (c j *2 . 0*pi*f r*pi* (4 . 0E-07) ) 

Etx= (377 . 0/ (2 . 0*pi) ) *cexp (c j *K2*d) *cos (Theta) *Ptx 
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Ety- (377.0/ (2.0*pi))*cexp(cj*K2*d)*cos(Theta)*Pty 
Ep x= (377 . 0/ (2 . 0*pi) ) *cexp ( c j *K2*d) *cos (Theta) *Ppx 
Epy= (377 . 0/ (2 . 0*pi) ) *cexp ( c j *K2*d) *cos (Theta) *Ppy 
ccccc 

Fl=dx*dy* (sin (Ky*dy2) / (Ky*dy2) ) * ( (sin (Kx*dx2) /(Kx*dx2) ) **2) 

F2=dx*dy* (sin (Kx*dx2) / (Kx*dx2) ) * ( (sin (Ky*dy2) / (Ky*dy2) ) **2) 

do 4 ip=l ,mm 

do 3 iq=l,nnpl 

irow= ( (ip-1) *nnpl) + iq 

Ph= cexp (c j * ( (-Kx*x ( ip) ) + (-Ky*y (iq) ) + (Ky*dy2) ) ) 

Vx(irow)= Di * ((Etx*Et)+(Epx*Ep)) *Fl*Ph 
if ((ip.eq.6) .and. (iq.eq.7)) Vmid=Vx(irow) 
ccccc write(6, 10)ip,iq,Vx(irow) 

3 continue 

4 continue 

do 7 ip=l,mmpl 

do 6 iq=l,nn 

irow= ( (ip-1) *nn) + iq 

Ph= cexp (c j * ( (-Kx*x ( ip) ) + (-Ky*y ( iq) ) + (Kx*dx2) ) ) 

Vy(irow)= Di * ((Ety*Et)+(Epy*Ep)) *F2*Ph 
ccccc write (6 , 10) ip , iq, Vy (irow) 

6 continue 

7 continue 

10 format (i3,i3,f 12.8, f 12.8) 

CCCCC 

CCCCC Now define shape function 
CCCCC 

Sx=l 

Sy=l 

do 550 m=i,mm 
do 540 n=l,nnpl 

dist=sqrt ( (x(m) **2)+( (y(n) -dy2) **2) ) 
if (dist.gt. 0.023) Sx(m,n)=0 
540 continue 
550 continue 

do 580 m=l ,mmpl 
do 570 n=l,nn 

dist=sqrt ( ( (x Cm) -dx2) **2) + (y (n) **2) ) 
if (dist.gt. 0.023) Sy(m,n)=0 
570 continue 
580 continue 

write(6,605) ((Sx(m,n) ,m=l,mm) ,n=nnpl , 1 ,-l) 

605 format (1211) 

write (6 , 606) ( (Sy (m,n) ,m=l ,mmpl) , n=nn , 1,-1) 

606 format (1311) 
irownew=0 

do 650 ip=l,mm 
do 640 iq=l,nnpl 
icolnew=0 

if (Sx(ip.iq) .eq.0) goto 640 
irow= ( (ip-1) *nnpl) + iq 
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ir ownew=irownew+ 1 

VN(irownew)=V(irow) 

do 620 m=l ,mm 

do 610 n=l,nnpl 

if (Sx(m,n) .eq.O) goto 610 

icol= ((m-l)*nnpl) + n 

icolnew=icolnew+l 

ZZNCirownew, icolnew)= ZZ(irow, icol) 

610 continue 
620 continue 

do 630 m=l,mmpl 

do 625 n=l ,nn 

if (Sy(m,n) . eq. 0) goto 625 

icol= C(m-l)*nn) + n + nmaxl 

icolnew=icolnew+l 

ZZNCirownew , icolnew) = ZZCirow , icol) 

625 continue 
630 continue 
640 continue 
650 continue 

do 700 ip=l,mmpl 
do 690 iq=l,nn 
icolnew=0 

if (Sy (ip, iq) . eq. 0) goto 690 

irow= ( (ip-1) *nn) + iq + nmaxl 

irownew=irownew+l 

VN (irownew)- V(irow) 

do 660 m= 1 ,mm 

do 655 n=l,nnpl 

if (Sx(m,n) . eq.O) goto 655 

icol= ((m-l)*nnpl) + n 

icolnew=icolnew+l 

ZZNCirownew, icolnew)= ZZ(irow, icol) 

655 continue 
660 continue 

do 680 m=l,mmpl 
do 670 n=l ,nn 
if (Sy (m,n) . eq.O) goto 670 
icol= ( (m-l)*nn) + n + nmaxl 
icolnew= icolnew+1 

ZZNCirownew, icolnew) = ZZCirow, icol) 

670 continue 
680 continue 
690 continue 
700 continue 
J=(0. 0,0.0) 

write (6 , 720) irownew , icolnew 
720 f ormat ( i4 , i4) 

750 call 12acg ( irownew, ZZN.nmax , VN, 1 , J ,f ac , ipvt ,wk) 
ccccc do 900 m=l,mm 

ccccc do 890 n=l,nnpl 
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ccccc if (Sx(m,n) .eq.O) goto 890 

ccccc mn=((m-l) *nnpl)+n 

ccccc write (6, 800)m, n, Jx(mn) , cabs (Jx(mn)) 

800 format (i3,i3,f 12.8, f 12.8, f 12.8) 

890 continue 
900 continue 


ccccc 

ccccc 

ccccc 

ccccc 

ccccc 


do 1000 m=l ,mmpl 
do 990 n=l,nn 
if (Sy(m,n) .eq.O) goto 990 
mn=( (m-1) *nn)+n 

write (6, 950)ra, n, Jy(mn) , cabs (Jy(mn)) 
950 format (i3,i3,f 12.8, f 12.8, f 12.8) 

990 continue 
1000 continue 
CCCCC 

CCCCC Now find RCS at angle (Theta, Phi) 

CCCCC 


irow=0 

Et=(0. 0,0.0) 

Ep=(0. 0,0.0) 

do 1100 ip=i,imn 
do 1050 iq=l,nnpl 
if (Sx(ip,iq) .eq.O) goto 1050 
irow= irow+1 

if ((ip*eq.6) .and. (iq.eq. 7)) Jmid=J(irow) 

Ph= cexp(cj * ( (-Kx*x(ip) ) + (-Ky*y (iq) ) + (Ky*dy2) ) ) 
Et= Et+(J(irow) *Etx*Fl*Ph) 

Ep= Ep+(J(irow)*Epx*Fl*Ph) 

1050 continue 
1100 continue 

do 1200 ip=l,mmpl 
do 1150 iq=l,nn 
if (Sy(ip,iq) .eq.O) goto 1150 
irow= irow+1 

Ph= cexp (c j * ( (-Kx*x ( ip) ) + (-Ky*y (iq) ) + (Kx*dx2) ) ) 
Et= Et+ ( J (irow) ^Ety*F2*Ph) 

Ep= Ep+(J(irow) *Epy^F2*Ph) 

1150 continue 
1200 continue 

1300 format (i3, i3,f 12.8, f 12.8) 

RCS= 4 . 0*pi*cabs (Et) *cabs (Et) 

RCS= 10. 0*alogl0(RCS) 

write (6 ,1500) (fr/(l .0E+09)) ,RCS, (Jmid/Vmid) 

1400 format (’Fr=’ ,f6.2, » RCS = \fl2.6,’ dBsm>) 

1500format(f8.4,f 12.6,2x,el2.6,2x,el2.6) 

2000 continue 
stop 
end 

ccccc 
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Solution for a Given Incidence Angle as a Function of Frequency 

The following program reads impedance matrix data files and computes backscatter at a 
given incidence angle as a function of frequency. Quadratic interpolation is used to find the 
impedance matrix if the frequency of interest is between that of t wo data files. The shape 
function is included to ensure that only the appropriate modes are retained in the solution. 

program BANDWIDTH 

parameter (mm=12 ,nn=12,mmpl=(mm+l) ,nnpl=(nn+l)) 

parameter (nmaxl=(mm*nnpl) ,nmax2=(mmpl*nn) ,nmax=(nmaxl+nmax2) ) 
parameter (NF=9) 

real Lx,dx,Ly, dy.Ko.fr ,x(mmpl) ,y(nnpl) 

real pi ,dx2 ,dy2 .Theta, Phi .Kx.Ky ,K2 ,F1 ,F2 

real f (NF) ,ZR(NF,nmax,nmax) ,ZI(NF,mnax,iimax) .epo.muo 

real tempr (nmax.nmax) .tempi (nmax.nmax) ,pxl ,px2,pyl ,py2 

real deltx, deity, ax, facty.factx.Tq.Tp.Tml ,Tn 

integer ipvt (nmax) .Sx(mm.nnpl) ,Sy(mmpl ,nn) 

complex er ,Ke ,Di ,Et ,Ep,Ph 

complex cj , K1 , Tm , Te , Ptx , Pty ,Ppx , Ppy ,Etx , Ety ,Epx ,Epy 
complex fac (nmax , nmax) .wk(nmax) .ass.bxx.res 
complex ZZ (nmax , nmax) , V (nmax) , RR (nmax , nmax) 
complex J (nmax) , JN 

complex ZZN( nmax, nmax) .VN(nmax) , Jmid.Vmid 
character* 12 fn 

DIMENSION U1 (3) ,U2(10) ,R1(3) ,R2(10) ,U(13) ,R(13) 

EQUIVALENCE (Ul(l) ,U(1)) , (U2(l) ,U(4)) , (Rl(l) ,R(1)) , (R2(l) ,R(4)) 

DATA Ul/. 11270166537925, .5, . 88729833462074/ ,U2/ . 01304673574141 , .06 

1746831665550 . . 16029521585048 . . 28330230293537 . . 42556283050918 . . 5744 
23716949081,-71669769706462, .83970478414951, .93253168334449, .986953 
326425858/ , Rl/ . 27777777777777 , . 44444444444444 , . 27777777777777/ , R2/ . 

403333567215434. . 07472567457529 . . 10954318125799. . 13463335965499 . . 14 

5776211235737. . 14776211235737 . . 13463335965499. . 10954318125799 .. 0747 

62567457529. . 03333567215434/ 

C 

do 20 IF=1,NF 
if (IF.eq.l) fn=’IMP6.0’ 
if (IF.eq.2) fn=’IMP6.5’ 
if (IF.eq.3) fn=’IMP7.0’ 
if (IF.eq.4) fn=’IMP7.5’ 
if (IF.eq.5) fn=’IMP8.0’ 
if (IF.eq.6) fn=’IMP8.5’ 
if (IF.eq.7) fn=’IMP9.0’ 
if (IF.eq.8) fn=’IMP9.5’ 
if (IF.eq.9) fn=’IMP10.0’ 
if (IF.eq.10) fn=’IMP10.5’ 
if (IF . eq. 11) fn=’IMP11.0’ 
if (IF . eq. 12) fn=’IMP11.5’ 
if (IF . eq. 13) fn=’IMP12.0’ 
if (IF . eq. 14) fn=’IMP12.5’ 
if (IF.eq. 15) fn=’IMP13.0’ 
if (IF.eq. 16) fn=’IMP13.5’ 
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if (IF. eq. 17) fn= , IMP14.0 J 

open ( 12 , f ile=fn , st atus= * OLD ’ , f orm= ’ UNFORMATTED ’ ) 

read ( 1 2 ) f r , Lx , dx , Ly , dy , d , Ko , Ke , er , x , y 

read(12)ZZ 

close (12) 

f (IF) = fr 

do 10 m=l , nmax 

do 10 n= 1 , nmax 

ZR(IF,m,n)= real (ZZ(m > n) ) 

ZI(IF,m,n)= aimag(ZZ(m,n) ) 

10 continue 
20 continue 
dx2=dx/2. 0 
dy2=dy/2 . 0 
pi= 2 . 0*asin(l . 0) 
cj= (0.0, 1.0) 
epo= 8.854E-12 
muo= pi*(4.0E-07) 

CCCCC 

CCCCC Start Freq. Loop 
CCCCC 

do 2000 IF=2000 , 6000 , 25 
f r= (float (IF)/ 1000 . 0)*(1 . 0E+09) 

Ko= 2 . 0*pi*f r*sqrt (epo*muo) 

Ke= 2 . 0*pi*f r*csqrt (er*epo*muo) 
ipt= -1 

call inter (NF , NF , f , nmax , ZR , 2 , f r , tempr , ipt , ierr ) 
ipt= -1 

call inter (NF,NF, f ,nmax, ZI , 2 , f r , tempi , ipt , ierr) 

ZZ= cmplx (tempr , tempi) 

CCCCC 

CCCCC Find Vx and Vy excitation for angle Theta and Phi 
CCCCC 

Theta= 60 . 0* (pi/180 . 0) 

Phi= 0 . 01*(pi/180 .0) 

Et= 1.0 
Ep= 0.0 

Kx= -Ko*sin(Theta)*cos(Phi) 

Ky= -Ko*sin(Theta) *sin(Phi) 

Kl= csqrt(cmplx(Ke**2)-(Kx**2)-(Ky**2) ) 
if (aimag(Kl) .gt .0.0) Kl= conjg(Kl) 

K2= Ko*cos (Theta) 

Tm= (er*K2*ccos(Kl*d) ) + (cj *Kl*csin(Kl*d) ) 

Te= (Kl*ccos (Kl*d) ) + (cj*K2*csin(Kl*d) ) 

Ptx= (Kl*Ko*cos (Phi) *csin(Kl*d) )/Tm 
Pty- (Kl*Ko*sin(Phi)*csin(Kl*d) )/Tm 
Ppx= - (Ko*Ko*sin(Phi) *csin (Kl*d) ) /Te 
Ppy= (Ko*Ko*cos (Phi) *csin(Kl*d) )/Te 
ccccc 

Di = -(4.0*pi)/(cj*2,0*pi*fr*muo) 

Etx= (377 . 0/ (2 . 0*pi) ) *cexp (c j *K2*d) *cos (Theta) *Ptx 
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Ety= (377 . 0/ (2 . 0*pi) ) *cexp (cj *K2*d) *cos (Theta) *Pty 
Epx= (377 . 0/ (2 . 0*pi) ) *cexp(c j *K2*d) *cos (Theta) *Ppx 
Epy= (377 . 0/ (2 . 0*pi) ) *cexp ( c j *K2*d) *cos (Theta) *Ppy 
ccccc 

Fl=dx*dy* (sin (Ky*dy2) / (Ky*dy2) ) * ( (sin (Kx*dx2) / (Kx*dx2) ) **2) 
F2=dx*dy* (sin (Kx*dx2) / (Kx*dx2) ) * ( (sin(Ky*dy2) / (Ky*dy2) ) **2) 
do 455 ip=l,mm 
do 454 iq=l,nnpi 
irow= ((ip-l)*nnpl) + iq 

Ph= cexp (c j * ( (-Kx*x (ip) ) + (~Ky*y (iq) ) + (Ky*dy2) ) ) 

V(irow)= Di * ( (Etx*Et) +(Epx*Ep) ) *Fl*Ph 
if ((ip.eq.6) .and. (iq.eq.7)) Vmid=V(irow) 
ccccc write (6,470) ip , iq,V(irow) 

454 continue 

455 continue 

do 465 ip=l,mmpl 
do 460 iq=l , nn 

irow= ((ip-l)*nn) + iq + nmaxl 

Ph= cexp (c j * ( (-Kx*x ( ip) ) + (-Ky*y ( iq) ) + (Kx*dx2) ) ) 

V(irow)= Di * ((Ety*Et)+(Epy*Ep)) *F2*Ph 
ccccc if ((ip.eq.6) .and. (iq.eq. 6)) Vmid=V(irow) 

ccccc write(6,470)ip,iq,V(irow) 

460 continue 
465 continue 

470 format (i3, i3,f 12 . 8,f 12 .8) 

CCCCC 

CCCCC Now define shape function 
CCCCC 
500 Sx=l 
Sy=l 

do 550 m=l,nmi 
do 540 n=l,nnpl 

dist=sqrt ( (x(m) **2) + ( (y (n)-dy2) **2) ) 
if (dist.gt. 0.023) Sx(m,n)=0 
540 continue 
550 continue 

do 580 m=l,mmpl 
do 570 n= 1 , nn 

dist=sqrt ( ( (x(m)-dx2)**2)+(y (n) **2)) 
if (dist.gt .0.023) Sy(m,n)=0 
570 continue 
580 continue 

if (IF. eq. 2000) write (6 , 605) ( (Sx(m,n) ,m=l ,mm) ,n=nnpl , 1 , -1) 

605 format (1211) 

if (IF. eq. 2000) write (6,606) ((Sy(m,n) ,m=l ,nunpl) ,n=nn, 1 ,-l) 

606 format (1311) 
irownew=0 

do 650 ip=l,mm 
do 640 iq=l ,nnpl 
icolnew=0 

if (Sx (ip , iq) . eq. 0) goto 640 
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irow= ((ip-l)*nnpl) + iq 
irownew=irownew+ 1 
VN (irownew) =V(irow) 
do 620 m=l , mm 
do 610 n=l,nnpl 
if (Sx(m,n) . eq. 0) goto 610 
icol= ((m-l)*nnpl) + n 
icolnew=icolnew+l 

ZZN (irownew, icolnew)= ZZ(irow, icol) 

610 continue 
620 continue 

do 630 m= 1 , mmp 1 

do 625 n=l,nn 

if (Sy(m,n) .eq. 0) goto 625 

icol= ((m-l)*nn) + n + nmaxl 

icolnew=icolnew+l 

ZZN (irownew, icolnew)= ZZ(irow, icol) 

625 continue 
630 continue 
640 continue 
650 continue 

do 700 ip=l , mmpl 
do 690 i q= 1 , nn 
icolnew=0 

if (Sy (ip, iq) . eq. 0) goto 690 
irow= ((ip-l)*nn) + iq + nmaxl 
irownew=irownew+ 1 
VN (irownew) = V(irow) 
do 660 m=l ,mm 
do 655 n=l,nnpl 
if (Sx(m,n) .eq.0) goto 655 
icol= ((m-l)*nnpl) + n 
icolnew=icolnew+l 

ZZN ( irownew, icolnew)= ZZ(irow, icol) 

655 continue 
660 continue 

do 680 m=l,mmpl 

do 670 n= 1 , nn 

if (Sy (m,n) . eq. 0) goto 670 

icol= ((m-l)*nn) + n + nmaxl 

icolnew= icolnew+1 

ZZN (irownew, icolnew) = ZZ(irow, icol) 

670 continue 
680 continue 
690 continue 
700 continue 
J=(0. 0,0.0) 

750 call 12acg (irownew , ZZN , nmax , VN , 1 , J , f ac , ipvt , wk) 
CCCCC 

CCCCC Now find RCS at angle (Theta, Phi) 

CCCCC 
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irow=0 

Et=(0. 0,0.0) 

Ep=(0. 0,0.0) 

do 1100 ip=l ,mm 

do 1050 iq=l,nnpl 

if (Sx(ip, iq) . eq. 0) goto 1050 

irow= irow+1 

ccccc write(6, 1040) f loat (ip) ,float(iq) , cabs (J(irow) ) 

1040 format (f8.2,f8.2,fl4.5) 

if ((ip.eq.6) .and. (iq.eq.7)) Jmid=J(irow) 

Ph= cexp(cj* ( (-Kx*x(ip) )+(-Ky*y(iq) )+(Ky*dy2) ) ) 

Et= Et+ ( J(irow) *Etx*Fl*Ph) 

Ep= Ep+( J(irow) *Epx*Fl*Ph) 

1050 continue 
1100 continue 

do 1200 ip=l,mmpl 
do 1150 iq=l,nn 
if (Sy (ip, iq) . eq. 0) goto 1150 
irow= irow+1 

ccccc write (6, 1040) float (ip) , float (iq) , cabs ( J(irow) ) 

ccccc if ((ip.eq.6) .and. (iq.eq. 6)) Jmid=J(irow) 

Ph= cexp(cj * ( (-Kx*x(ip) )+(-Ky*y (iq) )+(Kx*dx2) ) ) 

Et= Et+ ( J(irow) *Ety*F2*Ph) 

Ep= Ep+ ( J(irow) *Epy*F2*Ph) 

1150 continue 
1200 continue 

1300 format (i3,i3,f 12. 8, f 12. 8) 

RCS= 4.0*pi*cabs(Et)*cabs(Et) 
ccccc RCS= 4 . 0*pi*cabs (Ep) *cabs (Ep) 

RCS= 10 . 0*alogl0(RCS) 

JN=Jmid/Vmid 

write (6 , 1500) (f r/ (1 . 0E+09) ) ,RCS, JN 
1400 f ormat ( * Fr= * , f 6 . 2 , ’ RCS - dBsm 1 ) 

1500 format (f8.4,fl2.6,2x,el2.6,2x,el2.6) 
ccccc write (6 , 1600) Vmid, Jmid 

1600 format (4el4 . 6) 

2000 continue 
stop 
end 
ccccc 
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(a) Patch on a grounded slab of infinite extent. 
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Figure 2. Comparison between subdomain, entire-domain, and measured scattering from a 
rectangular microstrip patch antenna. L x = 1.88 cm; L y = 1.30 cm; d = 0.158 cm; e r — 2.17; 

Loss tangent = 0.001; 0\ ^ = 60° , 45°. 
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Figure 3. Comparison between subdomain and entire-domain calculated scattering from a 
circular microstrip patch antenna. Patch radius = 2.30 cm; d = 0.159 cm; £> = 2.20; 
Loss tangent = 0.0009; 6\ ft = 60°, 0°. 
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Figure 4. Comparison between subdomain; entire-domain, and measured scattering from a 
circular microstrip patch antenna. Patch radius = 0.71 cm; d — 0.07874 cm; e r = 2.20; 
Loss tangent = 0.0009; 9 l , S l = 63°, 0°. 
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Figure 5. Comparison between subdomain and measured scattering from an equilateral triangle 
microstrip patch antenna. Triangle side = 1.4 cm; d = 0.07874 cm; £ r = 2.33; Loss tangent 

= 0.001; e\ 4> l =60°, 180°. 
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Figure 6. Comparison between subdomain and measured scattering from a trapezoidal 
microstrip patch antenna. L x = 0.7 cm; L y = 1.4 cm; d = 0.07874 cm; e r = 2.33; Loss 
tangent = 0.001; 6 l , 4> l = 60°, 180°. 
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